// snow_qhc 合成ims气候场文件
// 像素值为10的时候是积雪


#include <string>
#include <stdio.h>
#include <string.h>
#include <iostream>
#include <fstream>
//#include "ogr_spatialref.h"//17
#include "gdal_priv.h"
#include <ctime>
#include <fstream>
#include "../sharedcodes-87/wftools.h"
#include "../sharedcodes-87/ajson5.h"
#include "../sharedcodes-87/wgt.h"


using namespace std;
using namespace ArduinoJson;
using std::cout;


//多期合成分类	IMS原始值
//范围外值	-9999
//缺测	-999
//常年陆地 -900 这个值仅在距平产品中有定义，用于区别距平为0的地区
//一次雪都没有且出现水的情况多于陆地和冰	-907
//一次雪都没有且出现陆地的情况多于水和冰	0（陆地颜色）
//一次雪都没有且出现冰的情况多于水和陆地	-909
//出现一次积雪及以上的像素	>0
#define MDAYS_PIXVAL_OUTSIDE -9999
#define MDAYS_PIXVAL_NODATA -999
#define MDAYS_PIXVAL_LAND_JP   -900
#define MDAYS_PIXVAL_WATER -907
#define MDAYS_PIXVAL_ICE -909













int main(int argc, char** argv)
{
	std::cout << "compute ims24km qhc 2021-6-7."<<endl ;
	std::cout<<"ims24km_qhc"<<endl ;


	//1997 -2020
	cout<<"gdal register"<<endl ;
	GDALAllRegister();

	string indir = "/home/hadoop/IdeaProjects/ims_asc2tif/tiff/" ;
	//string outfile = "ims24km_qhc00000301_00000531.tif";
	int startmmdd = 301;
	int endmmdd =   531;
	// int nyear = 2020 - 1997 +1 ;
	// int startyear = 1997 ;
	// int endyear = 2020  ;

	string outfile = "ims24km_20200301_20200531.tif";
	int nyear = 1 ;
	int startyear = 2020 ;
	int endyear = 2020  ;

	vector<string> allfiles ;
	//
	bool dirok = wft_direxist( indir ) ;
	if( dirok )
	{
		wft_get_allSelectedFiles( indir , "ims", ".tif", -1, "", allfiles);
	}else{
		cout<<"Error : indir is not exists."<<indir<<endl ;
		return 12;
	}

	int n = allfiles.size();


	wft_descSortStrVec( allfiles) ;


	vector<int> filldays , seadays, landdays , icedays , snowdays ;
	char proj[2048] ;
	double trans[6] ;
	bool projok =false ;
	int wid = 0 ;
	int hei = 0 ;
	int asize = 0 ;

	for (int i = 0; i < n; ++i )
	{//逐个文件处理
		string filepath = allfiles[i];
		string filename = wft_base_name(filepath);

		string ymd7 = filename.substr(3,7) ;
		int ymd7i = atof(ymd7.c_str()) ;
		int ymd8 = wft_convertymd7toymd8(ymd7i) ;

		cout<<ymd7<<"->"<<ymd8<<endl ;

		int year1 = ymd8 / 10000 ; 
		int mmdd =  ymd8 % 10000 ; 


		if( year1>=startyear && year1<=endyear && mmdd >= startmmdd && mmdd <= endmmdd ) {
			cout<<filename<<endl ;
			
			GDALDataset* inds1 = (GDALDataset*) GDALOpen(filepath.c_str() , GA_ReadOnly ) ;
			if( inds1 != 0 ){
				wid = inds1->GetRasterXSize() ;
				hei = inds1->GetRasterYSize() ;
				if( projok==false ){
					inds1->GetGeoTransform(trans);
					strcpy( proj ,  inds1->GetProjectionRef() ) ;
					asize =wid*hei  ;
					filldays.resize(asize,0) ;
					seadays.resize(asize,0);
					landdays.resize(asize,0);
					icedays.resize(asize,0);
					snowdays.resize(asize,0);

					projok = true ;
				}

				vector<short> tempindata(asize) ;
				//读取当前文件的数据
				inds1->GetRasterBand(1)->RasterIO(GF_Read, 0, 0, wid, hei,
					 tempindata.data() , wid, hei, GDT_Int16, 0, 0, 0);
				for(int it = 0 ; it < asize ; ++ it )
				{
					int iv = tempindata[it] ;
					if( iv == 0 ) {
						filldays[it] ++ ;
					}else if( iv == 1 ){
						seadays[it] ++ ;
					}else if( iv ==2 ){
						landdays[it]++ ;
					}else if( iv==3 ){
						icedays[it] ++ ;
					}else if( iv == 4 ){
						snowdays[it] ++ ;
					}
				}
				GDALClose(inds1) ;

			}else{
				cout<<"bad"<<endl ;
			}

		}
	}

	//average
	for(int it = 0 ; it < asize ; ++ it )
	{
		if( snowdays[it] > 0 ){
			snowdays[it] = snowdays[it] * 1.f / nyear + 0.5 ;
		}else {
			if( filldays[it] > seadays[it] && filldays[it] > landdays[it] && filldays[it] > icedays[it] ){
				snowdays[it] = -10 ;
			}else if( seadays[it] > filldays[it] && seadays[it] > landdays[it] && seadays[it] > icedays[it] )
			{
				snowdays[it] = -11 ;
			}else if( landdays[it] > filldays[it] && landdays[it] > seadays[it] && landdays[it] > icedays[it] )
			{
				snowdays[it] = -12 ;
			}else{
				snowdays[it] = -13 ;
			}

		}
	}



	cout<<"saving  qhc ..."<<endl ;

	GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("GTiff");
	//temp buffer

	GDALDataset* outds = driver->Create( outfile.c_str() , wid, hei, 1, GDT_Int16, 0);
	outds->SetGeoTransform(trans);
	outds->SetProjection(proj);
	outds->GetRasterBand(1)->RasterIO(GF_Write, 0, 0,wid, hei,
		 snowdays.data() , wid, hei, GDT_Int32,  0, 0, 0);			
	GDALClose(outds) ;


	std::cout << "All done." << std::endl;

	return 0;
}

